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Abstract 

The  homogenized  energy  model  (HEM)  is  a  unified  framework  for  modeling  hysteresis  in  ferroelectric, 
ferromagnetic,  and  ferroelastic  materials.  The  HEM  framework  combines  energy  analysis  at  the  lattice 
level  with  stochastic  homogenization  techniques,  based  on  the  assumption  that  quantities  such  as  inter¬ 
action  and  coercive  fields  are  manifestations  of  underlying  densities,  to  construct  macroscopic  material 
models.  In  this  paper,  we  focus  on  the  homogenized  energy  model  for  shape  memory  alloys  (SMA). 
Specifically,  we  develop  techniques  for  estimating  model  parameters  based  on  attributes  of  measured 
data.  Both  the  local  (mesoscopic)  and  macroscopic  models  are  described,  and  the  model  parameters’ 
relationship  to  the  material’s  response  are  discussed.  Using  these  relationships,  techniques  for  estimating 
model  parameters  are  presented.  The  techniques  are  applied  to  constant-temperature  stress-strain  and 
resistance-strain  data.  These  estimates  are  used  in  two  manners.  In  one  method,  the  estimates  are 
considered  fixed  and  only  the  HEM  density  functions  are  optimized.  For  SMA,  the  HEM  incorporates 
densities  for  the  interaction  and  relative  stress,  the  width  of  the  hysteresis  loop.  In  the  second  method, 
the  estimates  are  included  in  the  optimization  algorithm.  Both  cases  are  compared  to  experimental  data 
at  various  temperatures,  and  the  optimized  model  parameters  are  compared  to  the  initial  estimates. 


1  Introduction 

Shape  memory  alloys  (SMA)  are  novel  materials  that  exhibit  two  useful  characteristics.  The  super-elastic 
effect  allows  the  material  to  undergo  significant  strains  without  experiencing  plastic  deformation.  The  shape 
memory  effect  makes  the  material  capable  of  recovering  large  strains  upon  heating,  allowing  for  its  use  as 
an  actuator.  Furthermore,  SMA’s  bio-compatibility  and  high  energy  density  make  the  material  suitable  for 
in-vivo  actuation  applications,  potentially  revolutionizing  biomedical  devices. 

SMA  have  been  successfully  incorporated  into  a  number  of  prototypes  and  commercial  applications. 
The  super-elastic  effect  has  been  utilized  in  orthodontic  wires,  eye-glass  frames,  stents,  and  annuloplasty 
bands  [23].  Applications  using  the  shape  memory  effect  include  robotic  catheters  [33,34],  robotic  hands  [22], 
jet  chevrons  [9,10,26],  and  smart  inhalers  [21].  The  design  and  control  of  these  prototypes  are  complicated  by 
the  material’s  nonlinear,  hysteretic  dependence  on  stress  and  temperature.  This  complex  behavior  arises  from 
the  presence  of  three  crystalline  structures  (or  phases)  in  uniaxial  applications.  Computationally  efficient, 
accurate  models  of  the  material’s  behavior  are  necessary  for  design  optimization  and  control  algorithms. 
Furthermore,  the  model  adoption  is  greatly  facilitated  by  easily  identifying  parameters,  reducing  the  time 
required  to  achieve  accurate  model  fits. 

The  homogenized  energy  model  (HEM)  is  a  unified  framework  for  modeling  SMA,  ferroelectric,  and 
ferromagnetic  materials.  The  energy  origin  of  the  model  was  originally  investigated  for  SMA  by  Muller 
and  Wilmanksi  [17]  with  further  work  by  Achenbach,  Seelecke  [24,25],  and  Huo  [14].  The  original  mod¬ 
els  determined  the  equilibrium  phase  using  the  Gibbs  energy  to  predict  the  mesoscopic  (or  single-crystal) 
behavior.  The  HEM  framework  for  SMA  builds  on  these  models  by  homogenizing  (or  averaging)  the  local 
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response  to  predict  the  macroscopic  behavior  of  polycrystalline  materials.  The  homogenization  is  accom¬ 
plished  by  assuming  that  certain  parameters  are  manifestations  of  underlying  densities.  In  this  manner,  the 
model  quantifies  the  inherent  polycrystalline  nature  of  SMA  and  other  smart  materials  and  is  suitable  for 
applications  requiring  the  macroscopic  actuator  response. 

The  application  of  the  HEM  framework  to  SMA  was  originally  reported  in  [15],  and  the  model  has  also 
been  applied  to  ferroelectric  [13,29,30,32]  and  ferromagnetic  [19,27,28]  materials.  A  complete  description  of 
the  unified  framework  is  presented  in  [26,31].  Additionally,  the  HEM  has  been  incorporated  into  a  variety  of 
control  applications  [19,20].  A  computationally  efficient  implementation  of  the  model  is  described  in  [2,3]. 

In  this  paper,  we  focus  on  data-driven  techniques  for  estimating  model  parameters.  The  relationship 
between  the  parameters  and  experimentally  obtained  stress,  strain,  and  resistance  data  at  various  constant 
temperatures  is  discussed.  These  initial  parameter  estimates  are  used  in  two  manners.  In  one  method, 
the  parameters  are  considered  fixed  while  the  density  functions  are  optimized.  In  the  second  method,  the 
parameters  are  incorporated  into  the  optimization  process  with  the  density  functions.  Similar  techniques  for 
estimating  parameters  for  ferroelectric  and  ferromagnetic  models  are  presented  in  [11,12]. 

The  remainder  of  this  paper  is  organized  as  follows.  The  SMA  model  is  summarized  in  Section  2, 
including  the  local  model  based  on  the  Gibbs  energy  and  the  homogenization  techniques  used  to  construct 
the  macroscopic  model.  Techniques  for  identifying  the  model  parameters  from  experimental  data  are  detailed 
in  Section  3.  The  two  optimization  problems  are  presented  in  Section  4.  The  experimental  setup  is  described 
in  Section  5,  and  experimentally  obtained  data  is  used  to  estimate  the  initial  model  parameters.  Finally,  the 
results  of  the  optimized  model  are  presented  and  compared  to  the  experimental  data  in  Section  6. 


2  Homogenized  Energy  Model  for  SMA 

The  complex  behavior  of  SMA  is  due  to  the  presence  of  three  crystalline  phases  in  uniaxial  applications: 
austenite  and  two  martensite  phases.  Austenite  ( A )  is  induced  by  high  temperatures,  and  the  martensitic 
phases  are  induced  by  high  stresses.  Martensite  plus  (M- b)  is  stable  under  tensile  stresses,  and  martensite 
minus  (M— )  is  stable  under  compressive  stresses.  The  interaction  of  these  phases  yields  the  material’s 
complex,  thermo- mechanical  response.  The  transition  diagram  for  SMA  is  shown  in  Figure  1.  Numerous 
models  have  been  proposed  to  quantify  the  behavior  of  SMA,  including  both  physical  and  phenomenological 
models.  An  overview  of  these  models  is  presented  in  [26]. 

2.1  Mesoscopic  Model 

The  homogenized  energy  model  uses  the  free  energy  model  to  derive  the  kernel  for  the  local  (mesoscopic) 
material  behavior.  The  free  energy  model  of  SMA  uses  phase  fractions  (xa,  %m+,  and  xm-)  to  quantify  the 
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Figure  1:  Uniaxial  SMA  phase  transitions. 
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proportion  of  the  material  in  each  of  the  three  phases.  The  equilibrium  phase  fraction  depends  on  the  Gibbs 
energy 


where  the  Helmholtz  energy 
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is  represented  using  piecewise  quadratic  polynomials.  The  five  polynomials  correspond  to  the  three  phases 
and  the  barriers  between  them.  In  (2),  EA  and  Em  are  the  elastic  moduli  of  austenite  and  martensite, 
respectively,  and  £t  is  the  maximum  recoverable  strain.  The  temperature-dependent  inflection  points  are 
given  by 


£a(T) 


*_a(T) 

Ea 


and 

,  m-aM  (T)+„ 

£m\-L  )  —  y - 1"  &T 

and  represent  the  regions  where  each  of  phase  is  stable,  as  shown  in  Figure  2. 
The  critical  transition  stresses  are  given  by 


&a{T)  =  crL  +  A cft  ( T  —  Tl) 


(3) 


and 

vm{T)  =  &a{T)  -  cfr, 

where  cfl  is  the  transition  stress  from  A  to  M+  at  temperature  Tl ,  Acr^  is  the  hysteresis  loop’s  temperature- 
dependence,  and  gr  is  the  relative  stress  difference  in  loading  and  unloading.  The  parameters  listed  in  Table  1 
can  be  related  to  experimental  data. 


Figure  2:  Relationship  between  the  inflection  points  for  austenite  and  martensite  and  the  critical  transition 
stresses. 
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Table  1:  SMA  model  parameters  in  the  Gibbs  energy  and  likelihood  relations. 


Variable 

Description 

Ea 

Austenite  elastic  modulus 

Em 

Martensite  elastic  modulus 

Martensite  transition  stress  at  Tl 

tl 

Lower  transition  temperature 

ctr 

Stress  hysteresis  loop  width 

Act  x 

Hysteresis  loop’s  temperature  dependence 

st 

Maximum  recoverable  strain 

T 

Relaxation  time 

V 

Layer  volume 

The  unknown  coefficients  in  (2)  are  found  by  ensuring  G 1  continuity,  yielding 

i/mx  _  EM  (st  -  £m(T))  +  Easa(T) 

[  2  (sa(T)  -  £m(T)) 

R/Tx  _  £a(T)  ( Ea  (_£t  -  £m(T))  +  EAeA(T)) 

[  1  eA{T)-eM{T) 

_  2A/3(T)  (eA(T)  -  eM(T))  +  £A(T)2  (Em  (eT  ~  £m(T ))  +  EA£m(T)) 
(  ’  2  (sa(T)  -  esM(T)) 

and 

Ap(T)  =  t  (Em  (st  ~  £a(T ))  (sT  ~  £m(T ))  -  EA£a(T)£M(T))  . 

The  evolution  of  the  phase  fractions  is  governed  by  the  coupled  differential  equations 

XM+(t)  =  ~p+AXM+(t)  +  PA+XA(t) 

XM-(t)  =  ~P-AXM-(t)  +pA-XA(t )• 

The  final  phase  fraction  is  found  from  the  conservation  equation 

1  =  XA(t)  +  XM+(t)  +  xM-(t). 


(4) 

(5) 


In  (4),  Pap  is  the  likelihood  of  transition  from  phase  a  to  ft  (a.  ft  =  M+,M—,  and  A).  The  transition 
rates  depend  on  the  Gibbs  energy  and  can  be  found  using  a  number  of  methods.  The  original  method 
developed  by  Achenbach  [1]  and  Muller  and  Seelecke  [16]  is  based  on  statistical  mechanics  and  the  sampling 
of  the  energy  barriers  by  the  equilibrium  phases.  Here,  we  assume  that  the  transition  rates  depend  on  the 
activation  energy,  or  the  barriers  in  the  energy  landscape.  Using  the  activation  energy  A the  transition 
rate  is 


Pa/3  (c,  T)  =  1  exp 


VAGaP  (a,  T) 
kT 


where  k  is  Boltzmann’s  constant. 

The  relaxation  time  r  and  layer  volume  V  must  be  identified  using  experimental  data  and  can  be  included 
in  the  optimization  of  model  parameters.  A  schematic  of  the  four  barriers  is  shown  in  Figure  3  for  a  small 
applied  stress  and  low  temperature.  The  barriers  A Gap  (cr,  T)  are  determined  by  finding  the  difference  in 
minima  and  maxima  of  the  polynomials  in  (1).  For  example, 


AG_ A  (a,  T)  =  G2  (a,  e*2 ,  T)  -  G1  (<r,  e^T) . 
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Figure  3:  Gibbs  energy  landscape  with  barriers  highlighted. 


where  and  s{  respectively  solve 
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The  strain  for  a  lattice  element, 

e(t)  =  xA(t)(eA)  +  xM+(t)(£M+)  +  xM-(t){eM-),  (6) 


depends  on  the  phase  fractions  and  their  respective  equilibrium  strains.  Assuming  low  thermal  activation, 
the  equilibrium  strains  are  given  by 


{sa)  = 


<r(t) 
Ea  ’ 


(eM+)  = 


Em 


(«-> = 


st- 


The  electrical  resistance  in  SMA  also  depends  on  the  phase  fractions.  In  some  applications,  for  example 
a  wire,  the  resistance  is  useful  for  self-sensing,  where  the  measured  resistance  is  used  to  predict  the  strain  in 
the  wire  [6-8].  As  detailed  in  [6],  the  resistivity 


Psma(T )  =  XA(t)pA(T )  +  [xM+(t)  +  XM-(T)]pM(T) 


depends  on  temperature,  where  the  resistivity  of  each  phase 


Pa{T)  —  poA  (1  +  aAAT) 

Pm(T)  =  Pom  (1  +  «mAT) 

includes  a  nominal  resistivity  pop  and  the  thermal  dependence  ap.  The  total  resistance  in  the  wire, 

C^f  r  _  Psma{T)L(s) 

1 l'£)~  A(ee)  ’ 

accounts  for  the  change  in  length  of  the  wire  due  to  strain 

L(s)  =  Lo  (1  +  s  (t)) . 


(7) 

(8) 

(9) 
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Additionally,  the  cross-sectional  area 


A(ee)  =  7T  (i?o  (1  -  V£e))2  , 


is  corrected  to  account  for  elastic  deformation,  where  the  elastic  strain  is 


£e  =  cr(t) 


f xM+(t)  ~  xM-(t) 
V  Em 


xA(t)\ 

Ea  )  ' 


Here,  To  is  the  austenitic  (zero  strain)  length  of  the  wire,  Ro  is  the  nominal  wire  radius,  and  v  is  Poisson’s 
ratio.  Similar  to  r  and  V,  the  resistance  model  parameters  (pom,  Poa ,  and  v)  must  be  estimated 

from  experimental  data  and  can  be  included  in  the  optimization  scheme  as  well. 

In  addition  to  the  phase  fraction  ODEs  (4)  and  the  constitutive  relations  (6)  and  (9),  an  internal  energy 
balance  is  necessary.  The  thermal  evolution  equation  at  the  local  level  is 


pcvT(t)  =  —h  (T(t)  -  T^)  -K(T(t)  -  Te(t ))  +j(t)  +  pH  (xM+(t)  +  xM-(t)) , 


where  p  is  the  SMA  density,  cv  is  the  specific  heat,  h  is  the  convective  heat  transfer  coefficient  with  ambient  air 
temperature  T^,  k  is  the  conduction  coefficient  with  surrounding  temperature  Te,  j(t)  is  the  electrical  input 
power,  and  H  is  the  latent  heat  associated  with  the  phase  transformations.  Here,  we  are  using  experimental 
data  collected  in  a  constant-temperature  water  bath  and  assume  that  the  temperature  change  due  to  the 
latent  heats  is  negligible  due  to  the  presence  of  a  large  thermal  mass  (water);  hence,  the  thermal  evolution 
equation  is  unnecessary. 

The  ODEs  for  the  phase  fractions  (4)  can  be  discretized  and  solved  using  an  implicit  Euler  scheme.  In 
the  discretized  form,  the  phase  fractions  xk+1  =  [x^;,  x^^]T,  at  time  step  k  +  1,  are  given  by 

xk+1  =  xk  +  A  t[A(tk+1)xkJrl  +  6(t/c+1)], 


where 

-P+A  (crk+1  ,Tk+1)  ~PA+  (ufc+1,Tfc+1)  -pA+  (ak+\Tk+1) 

-PA-  (crk+1,Tk+1)  -p-A  (ak+1,Tk+1)  ~  PA-  (crk+1,Tk+1) 

'pA+  (afc+\  Tk+l 
PA-  (ak+1,Tk+1)\ 

The  initial  condition  for  the  phase  fractions  depends  on  the  temperature  and  stress;  however,  initializing 
experiments  and  simulations  under  zero  stress  simplifies  the  initial  condition.  At  temperatures  above  the 
austenitic  start  temperature  of  SMA,  Tas  ,  the  initial  phase  fraction  is  assumed  to  be  fully  austenitic  and 
xm+( 0)  =  xm-{ 0)  =  0.  The  austenitic  start  temperature  is  the  approximate  temperature  at  which  SMA 
would  transition  to  austenite  under  zero  stress  and  is  available  from  manufacturer  data.  The  value  can 
also  be  determined  experimentally  by  straining  a  wire  and  observing  the  temperature  at  which  it  begins  to 
contract  under  zero  stress.  At  lower  temperatures,  experimental  results  have  shown  a  propensity  to  exist 
in  the  M-f  phase,  leading  to  a  non-zero  initial  strain  even  after  heating  and  cooling  under  zero  stress  [5]. 
This  effect  may  be  due  to  material  inhomogeneities  or  the  two-way  shape  memory  effect.  To  quantify  the 
non-zero  initial  strain,  the  initial  condition  is  based  on  this  strain  value.  At  low  temperatures,  the  initial 
strain  is  given  by 

xm+(0)^t  —  xm-(0)^t  =  e(0). 

Substituting  in  the  conservation  relation  (5)  and  setting  the  xa  fraction  to  zero  yields 

XM+{ o)  =  +  —  ,  XM-(0)  =  1  -  xM+(0). 

1st 


A  (tk+1)  = 
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2.2  Homogenized  Energy  Model 

The  mesoscopic  or  lattice-level  model  quantifies  the  single  crystal  behavior  of  SMA  by  assuming  uniform 
properties  and  disregarding  material  and  stress  inhomogeneities;  however,  SMA  is  inherently  polycrystalline. 
The  HEM  quantifies  the  macroscopic  behavior  by  assuming  that  the  relative  stress  a r  and  interaction  stress 
a i  are  manifestations  of  underlying  densities  vr{or)  and  ^7(07).  The  macroscopic  strain  includes  these 
effects  and  is  given  by 

POO  POO 

e(a(t),T(t))=  /  pR(aR)pI(aI)£(a(t) +aI,T(t),aR)daIdaR.  (10) 

J  0  J  —  oo 


The  relative  stress  density  changes  the  hysteresis  loop  width  (the  difference  in  the  loading  and  unloading 
stress,  as  shown  in  Figure  4(a)).  However,  all  kernels  switch  from  austenite  to  martensite  at  the  same  stress 
(point  A).  The  relative  stress  also  affects  the  height  of  the  barrier  between  M+  or  M—  and  A  in  the  Gibbs 
energy  landscape;  see  Figure  4(c). 

The  interaction  stress  density  affects  both  loading  (switching  from  austenite  to  martensite)  and  unloading 
(switching  from  martensite  to  austenite).  For  a  given  temperature  and  applied  stress,  the  interaction  stress 
density  shifts  the  hysteresis  loops  up  and  down  (Figure  4(b)),  but  all  kernels  have  the  same  width  and  the 
same  parameters.  As  illustrated  in  Figure  4(d),  this  behavior  corresponds  to  a  tilt  in  the  Gibbs  energy  since 
the  interaction  stress  changes  the  effective  stress. 

Both  densities  are  necessary  to  accurately  model  the  major  loop  and  minor  loop  behavior.  Without  the 
interaction  density,  all  kernels  switch  from  austenite  to  martensite  at  the  same  stress  (between  points  A 
and  B  in  Figure  5(a)).  The  only  changes  occur  on  unloading  as  the  material  switches  from  martensite  to 
austenite  (between  points  C  and  D  in  Figure  5(a))  over  a  range  of  stresses  corresponding  to  the  different 
relative  stresses. 

Without  the  relative  stress  density,  the  minor  loops  will  be  anhysteretic,  as  illustrated  in  Figure  5(b). 
The  kernels  begin  switching  from  austenite  to  martensite  at  point  A.  If  the  stress  is  decreased  from  point 
B  to  C,  all  of  the  kernels  that  switched  to  martensite  will  remain  in  martensite.  Since  all  kernels  have  the 
same  relative  stress  or ,  the  stress  will  have  to  decrease  to  approximately  point  D  before  switching  back  to 
austenite  will  occur.  Therefore,  a  minor  loop  between  points  B  and  C  will  be  anhysteretic  (hysteresis- free) , 
which  is  contrary  to  the  observed  behavior. 

The  relative  stress  and  interaction  stress  densities  can  take  multiple  forms,  subject  to  certain  conditions. 
Both  densities  must  be  integrable.  The  relative  stress  density  is  defined  only  for  or  >  0,  since  the  switch  from 
austenite  to  martensite  occurs  at  a  higher  stress  than  the  reverse  transition.  Furthermore,  the  interaction 
stress  density  is  assumed  to  be  symmetric.  In  [4],  we  compare  two  methods  for  representing  the  densities. 
One  method  uses  a  log-normal  probability  density  function  for  i 'r(or)  and  a  normal  density  for  vi(oi).  The 
other  method  uses  a  linear  combination  of  underlying  log-normal  and  normal  probability  density  functions. 
Here,  we  use  the  second  method  as  it  provides  greater  flexibility  and  accuracy.  The  two  densities  are  then 
expressed  as 


vr(or)  = 


Ma  Ka 

TEE  ,rn  4>k  (os) 
1  m=  1  k= 1 


Mo,  Ka 

Ci  =  EEa^ 

m=  1  k= 1 


(ii) 


and 


vi(oi)  = 


1 

G2 


(07) 

k=  1 


C2  =  ym 


k= 1 


(12) 


where  the  coefficients  ak:m  and  fik  are  identified  by  optimizing  the  model  response  to  observed  experimental 
data.  The  coefficients  C\  and  C2  ensure  that  the  densities  integrate  to  unity.  The  basis  functions  (\>k,m  (&r) 
and  vpk  (07)  are  given  by 


(0\r) 


- -=  exp  (-[In  (aR)  -  fj,Rj2/2c2k) 

ckaRV2Tr 


(13) 
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and 


1 


(14) 


ipk  (o'/)  = 


bkV^TT 


exp  {-a]/2bl). 


Initial  estimates  for  the  relative  stress  log-normal  mean  /hr  and  standard  deviation  c  and  interaction  stress 
standard  deviation  6  are  obtained  from  experimental  data.  The  initial  estimates  are  then  perturbed  to 
obtain  the  desired  basis  function  parameters  /iRm,  C&  and  6^.  For  example,  after  estimating  the  standard 
deviation  for  the  interaction  stress,  6,  four  additional  basis  function  are  created  to  have  standard  deviations  of 
0.256,  0.56,  26,  and  46.  The  parameters  fiRrn  and  c/~  are  perturbed  in  a  similar  manner  about  experimentally- 
determined  value  of  fiR  and  c. 

The  macroscopic  model  (10)  is  integrated  using  4-point  Gaussian  quadrature  on  20  equal  intervals.  The 
discretized  model  is  evaluated  using 


Ni  Nj 

£  ( <r(t),T(t ))  =  ^ ^  vr  ( aRi )  vi  (cr (t)  +  a/.,T(t),aRi )  wtWj, 


i= 1  3  =  1 


(15) 


where  the  weights  Wi  and  Wj  are  determined  by  Gaussian  quadrature.  The  summations  in  (15)  are  easily 


(a) 


Strain 


Strain 


(c)  (d) 

Figure  4:  Homogenenized  energy  model  kernels  for  different  values  of  (a)  cfr  and  (b)  aj.  Gibb’s  energy 
landscape  for  different  values  of  (c)  gr  and  (d)  cr/. 
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Figure  5:  Hysteresis  behavior  when  (a)  vr{gr)  is  distributed  but  i/j(crj)  is  constant,  and  (b)  z/j(crj)  is 
distributed  but  vr(cfr)  is  constant. 


represented  as  matrix  multiplications,  increasing  the  computational  efficiency  of  the  model.  In  matrix  form, 
the  discrete  model  is 

e(a(t),T(t))  =  VTTW. , 

where 

VT  =  —  [wii/r(<tRi),  ■  ■  ■  ,wNivR{aRN,)} 

«i 

and  i 

WT  =  ,W1 v^/(ct/n.)] 

incorporate  the  quadrature  weights  and  densities  evaluated  at  the  quadrature  points.  The  matrix 

[T]»i  =e{a(t)+(jIj,T(t),(TRi) 

is  the  strain  kernel  evaluated  at  the  quadrature  points. 

3  Data-driven  Techniques  to  Determine  Initial  Parameter  Esti¬ 
mates 

3.1  Mesoscopic  Model  Parameters 

One  of  the  advantages  of  the  homogenized  energy  model  is  that  the  model  parameters  listed  in  Table  1  can 
be  related  directly  to  experimental  data.  Stress,  strain,  resistance,  and  temperature  data  can  be  used  to 
estimate  these  parameters.  These  estimates  can  either  be  used  directly  as  the  model  parameters  or  as  initial 
estimates  for  optimization  algorithms.  Accurate  initial  parameter  estimates  greatly  enhance  optimization 
convergence.  Furthermore,  accurate  initial  conditions  may  help  overcome  the  likelihood  of  converging  to 
local  minima  instead  of  the  global  minimum. 

Many  of  the  parameters  can  be  directly  estimated  from  stress-strain  or  resistance-strain  data.  For  ex¬ 
ample,  the  elastic  moduli  Ea  and  Em  are  determined  from  high  temperature  (above  the  austenite  start 
temperature  Tas)  stress-strain  data,  as  shown  in  Figure  6.  The  austenitic  elastic  modulus  Ea  is  the  slope 
of  the  stress-strain  curve  as  the  SMA  wire  is  initially  loaded  (line  ab  in  Figure  6)  .  The  martensitic  elastic 
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Figure  6:  Relationship  between  the  stress-strain  response  of  SMA  and  model  parameters  Ea  (line  ab ),  Em 
(line  cd),  and  Et  (point  e). 


modulus  is  the  slope  of  the  stress-strain  curve  at  high  stresses  during  unloading  (line  cd  in  Figure  6).  The 
maximum  recoverable  strain  St  is  the  intersection  of  cd  with  the  zero  stress  axis  (point  e  in  Figure  6). 

The  stress-strain  hysteresis  loop  is  temperature-dependent  and  satisfies  the  relation 

0A  =  vl  +  Act t  ( T(t )  -  Tl)  , 

where  a  a  is  the  transition  stress  from  austenite  to  martensite.  The  transition  stress  cfl  and  Tl  are  coupled 
and  not  unique.  Therefore,  the  values  can  be  determined  from  the  stress-strain  response  above  Ta5,  as 
depicted  in  Figure  7.  At  a  given  temperature  TL  =  crL  is  the  stress  at  the  midpoint  of  the  transition 
from  austenite  to  martensite  (point  /  in  Figure  7).  The  transition  stress  from  austenite  to  martensite  varies 
linearly  with  temperature  based  on  A ctt-  An  initial  estimate  for  Act t  is  determined  from  the  stress-strain 
response  at  two  different  temperatures.  Given  the  stress  at  points  /  and  g  in  Figure  7,  it  follows  that 

Act  =  •  (16) 

The  relaxation  time  r  and  volume  layer  V  are  determined  from  creep  data.  Creep  data  is  obtained  by 
quickly  loading  the  SMA  to  a  specified  strain  £\  and  then  holding  the  strain  constant  (point  a  to  point  b  in 
Figure  8(a)).  While  the  strain  is  held  constant  at  e i,  the  stress  decreases  (Figure  8(c))  as  austenite  converts 
to  M+;  see  Figure  8(d).  Alternatively,  the  SMA  can  be  loaded  to  a  specified  stress  (a)  while  the  strain  is 
monitored. 

Assuming  that  during  relaxation,  the  phase  fraction  transitions  from  A  to  M+  only,  the  phase  fraction 
ODE  (4)  reduces  to 

&M+{t)  =  Pa+  (<r,T)  xA(t)  =  Pa+  (c,  T)  (1  —  xM+(t)) .  (17) 

Neglecting  interaction  and  material  inhomogeneity  effects  and  assuming  to  =  0  and  XM+(to)  =  xm+0 ,  the 
analytical  solution  to  (17)  is 


XM+(t)  =  {xM+o  +  1)  exp  ( PA+  (<7,  T )  t)  -  1. 


(18) 


In  (18),  the  transition  rate  pa+  (<t,T)  depends  on  the  the  stress;  however,  it  is  assumed  that  at  strains  near 
the  maximum  recoverable  strain  et  the  energy  barrier  A Ga+  (c,  T )  has  been  eliminated,  yielding 


pA+  ( cf ,  T)  =  -  exp 

T 


VAGa+  (<r,r) 
kbT 


1 

r 


(19) 
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The  assumption  that  the  energy  barrier  is  zero  eliminates  the  stress  dependence  and  the  volume  layer  V, 
leaving  the  unknown  model  parameter  r.  As  shown  in  Figure  9,  as  the  material  begins  switching  from  A  to 
M-h  (point  a),  the  energy  barrier  is  large  and  its  effect  on  the  transition  rate  cannot  be  neglected.  However, 
at  strains  near  St  (point  6),  the  energy  barrier  is  nearly  eliminated.  Furthermore,  we  note  that 


V 

kbT  <<  5 


providing  further  justification  for  the  assumption 


VAGa+(ct,T)\ 
kT  ) 


1 


at  strains  near  Et-  Rate  effects  in  this  region  are  attributed  solely  to  the  relaxation  time  r. 

Using  the  measured  creep  data,  we  assume  that  the  decay  in  stress  is  proportional  to  the  decay  in 
austenite  (or  the  increase  in  M+  since  the  ODEs  are  developed  using  M+  and  M— ).  Letting  d\  and 
denote  the  change  in  stress  shown  in  Figure  10(a)  and  using  the  assumption  that  the  decay  rate  for  stress  is 
equal  to  the  increase  in  M+,  we  obtain 

%M+{t 2)  ~  ffM+fe)  _  £3  ~  0~2  _  (20) 

XM+itl)*- XM+(t  2)  (J2—Cri  d\ 


Assuming  constant  time  steps  At  and  using  (18),  Equation  (20)  simplifies  to 

exp  (-2  •  At  •  pa+  0?  T))  -  exp  (-3  •  At  •  pA+  T))  _  d^ 
exp  (-At  •  pa+  (cr,  T))  -  exp  (-2  •  At  •  pA+  (cr,  T ))  d{ 


Solving  for  pA+  (&,  T)  yields 

Pa+  (ct,T)  =  ^  In 

Substituting  the  simplified  transition  likelihood  (19)  gives 


(21) 


(22) 


Figure  7:  Relationship  between  the  stress-strain  response  of  SMA  and  model  parameters  cjl  (point  /),  Tl 
(temperature  corresponding  to  stress  at  point  /),  and  Aar  (change  in  stress  observed  between  /  and  g  over 
change  in  temperature  X2  -  T\). 
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(a)  (b) 


(c)  (d) 

Figure  8:  Creep  data  used  to  obtain  estimates  of  r  and  V:  (a)  strain  input  (constant  from  a  to  6);  (b) 
resulting  stress-strain  response;  (c)  stress  decay  from  a  to  b  as  strain  is  held  constant;  (d)  assumed  increase 
in  M+  phase  fraction  as  strain  is  held  constant. 


After  solving  for  r,  creep  data  at  a  lower  strain  (point  a  in  Figure  9(a))  is  used  to  determine  the  volume 
layer  V.  The  same  process  is  used,  except 


Pa+  (cr,  T)  =  -  exp 

T 


(  V  AG  a+  (&,T)\  1  ,  {d,\ 

l — hf — J  =  iilnUJ 


(23) 


is  solved  for  V.  The  barrier  A Ga+  (cr,  T)  depends  on  the  stress;  however,  (23)  is  derived  assuming  the  stress 
is  constant  over  this  interval,  which  is  not  the  case.  The  stress  at  the  midpoint  can  be  used  or  Equation  (21) 
can  be  solved  for  V  using  the  stresses  at  £i,  £2,  and  £3. 

The  resistance  model  parameters  for  martensite  are  estimated  using  the  measured  resistance  at  different 
low  temperatures.  At  low  temperatures,  the  equilibrium  phase  is  M+  after  loading  and  unloading.  At  this 
point  (a  strain  of  et  and  point  a  in  Figure  11(a)),  effects  of  olm  and  v  are  negated.  The  resistivity  is 


Pm  0 


VtaltRl 

L{£e  J  ’ 


(24) 
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where  Qa  is  the  measured  resistance  and  Ro  is  the  SMA  radius.  The  Poisson’s  ratio  is  calculated  using  the 
low-temperature  resistance  at  two  different  strains  and  by  solving 


tA  —  =  pMO 


L{eeb) 


Tt(R0-  RoV£b) 


L(£ea)\ 

nRl  ) 


(25) 


for  v.  The  austenite  nominal  resistivity  is  determined  using  high-temperature  resistance  strain  data  shown 


Figure  9:  (a)  High-temperature  response  of  SMA  during  loading  using  r  =  0.1  s  and  V  =  5  x  10  23  m3;  (b) 
energy  barrier  between  A  and  M+  during  loading. 


Time  Time 

(a)  (b) 

Figure  10:  Estimation  of  r:  (a)  stress  decay  during  creep  and  change  in  stress  d\  and  d 2  over  equal  time 
intervals;  (b)  corresponding  increase  in  M+  phase  fraction  during  creep. 
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(a)  (b) 


Figure  11:  Estimation  of  resistance  model  parameters:  (a)  pMo  is  estimated  from  the  resistance  at  point 
a  and  v  estimated  using  the  change  in  resistance  between  points  b  and  a;  (b)  pao  using  high-temperature 
resistance-strain  data  (point  c). 


in  Figure  11.  Using  the  resistance  at  zero  strain,  the  resistivity  is 


Pao 


QcttRq 

We  J' 


(26) 


Note  that  this  calculation  assumes  that  the  temperature  dependence  is  based  on  a  different  nominal  tem¬ 
perature  than  the  martensite  dependence. 

As  shown  in  Figure  12,  the  temperature  dependence  of  the  resistivities  olm  and  a  a  can  be  found  by 
varying  the  temperature  under  zero  stress.  After  loading  and  unloading  at  low  temperatures,  the  wire  will 
be  in  M-f.  As  the  temperature  rises,  the  resistance  (and  resistivity)  rises  due  to  olm •  Near  the  austenite  start 
temperature,  the  resistance  drops,  due  to  the  lower  resistivity  of  austenite.  Further  increases  in  temperature 
lead  to  increases  in  resistance  due  to  a  a-  Using  the  points  indicated  in  Figure  12(a),  we  obtain 


OLM 


n2A(ee)/L(se)  -  n1A(ge)/£(ge) 
T2  -  Ti 


(27) 


and 


OtA 


Q,±A{ee)/L{ee)  -  fl3A(£e)/L(£e) 

T4-T3 


(28) 


Alternatively,  olm  and  a. a  can  be  identified  using  experimental  resistance-strain  data  at  two  different 
temperatures  after  identifying  pmo,  Pao ,  and  v.  For  example,  the  measured  resistance  at  d  and  e  in  Fig¬ 
ure  12(b)  can  be  used  to  find  a  a  using  (7).  Similarly,  the  measured  resistance  at  /  and  g  can  be  used  to 
find  olm  using  (8). 


3.2  Homogenized  Energy  Model  Parameters 

An  initial  estimate  for  the  interaction  stress  standard  deviation  b  is  obtained  from  high-temperature  stress- 
strain  data.  As  shown  in  Figure  5(b),  the  interaction  stress  causes  the  material  to  switch  from  A  to  M+ 
over  a  range  of  stresses  instead  of  at  a  a-  Using  this  range  of  stress  and  assuming  a  normal  distribution, 
approximately  95%  of  the  kernels  switch  from  A  to  M+  between  oa  —  2  b  and  oa  +  2  b.  Note  that  a  a  is 
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determined  from  (3)  after  estimating  cr^,  ?l,  and  Actt-  Switching  from  A  to  M+  begins  where  the  stress- 
strain  response  deviates  from  elastic  loading  along  Ea  and  ends  where  the  material  elastically  deforms  along 
Em]  see  Figure  13.  After  estimating  the  interaction  stress  standard  deviation,  the  interaction  stress  density 
kernels  (14)  are  initialized  with  multiples  of  6,  yielding  the  desired  number  of  b^. 

Using  the  switching  region  identified  in  Figure  13,  the  parameters  for  the  relative  stress  density  vr(<jr) 
can  be  estimated  using  characteristics  of  log-normal  densities  [18].  The  hysteresis  loop  width  gr  is  estimated 
using  the  difference  in  loading  and  unloading  at  the  midpoint  of  the  major  loop  (a a  on  loading),  as  shown 
Figure  14.  Assuming  this  value  is  the  mode,  the  log-normal  density  mean  is  given  by 

Hr  =  In  (aR) .  (29) 


(a)  (b) 

Figure  12:  Estimation  of  ctM  and  a  a  using  resistance-temperature  data:  (a)  resistance-temperature  data  at 
zero  stress;  (b)  resistance-strain  data  at  two  different  temperatures  (points  d  and  e  for  a  a  and  points  /  and 
9  for  aM)- 


Strain 


T2>T! 


Figure  13:  Relationship  between  the  stress-strain  response  of  SMA  and  the  initial  estimate  for  the  interaction 
stress  standard  deviation  b. 
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On  unloading,  switching  ends  when  z/j(<7j)  stops  intersecting  with  vr(gr)  at  the  95.4%  confidence  interval. 
Using  this  assumption  and  identifying  the  end  of  switching  ( gr  in  Figure  14),  the  standard  deviation  is 


c 


1 

2 


(In  ( gf  -  2b)  -  fiR)  . 


(30) 


The  relative  stress  density  basis  functions  (13)  are  initialized  with  multiples  of  the  macroscopic  estimate  for 
the  log-normal  mean  /ir  and  standard  deviation  c,  which  yields  the  desired  values  for  /jjrtu  and  c&. 

An  alternative  method  to  identify  the  density  parameters  is  to  fit  a  log-normal  and  normal  density  directly 
to  the  experimental  data.  After  identifying  a  a  and  the  switching  region  in  Figure  13,  the  interaction  stress 
standard  deviation  is 


b  = 


Nk 


\  Nk 


-  y]  ((Tfc  -  aAf, 


(31) 


where  crk  is  the  experimentally  measured  stress  within  the  switching  region. 
The  experimentally  observed  relative  stress 


-  k  _  Ak 
aR  —  aA 


is  estimated  at  each  measured  strain  ek  by  taking  the  difference  in  the  loading  stress  crkA  and  unloading 
stress  .  This  method  determines  the  experimentally  observed  relative  stress  throughout  the  entire  major 
loop,  as  shown  in  Figure  15.  Using  this  collection  of  estimates,  the  mean  E  (<r|Q  and  variance  Var  (b^)  are 
calculated.  Again  using  the  characteristics  of  log-normal  distributions,  the  density  parameters  are  given  by 


and 


Hr  =  In  (E(4))-^ln 


c2  =  In 


Varp|)\ 

e(4)2;' 


(32) 


(33) 


The  methods  for  estimating  the  model  parameters  are  summarized  in  Table  2.  The  methods  described 
in  this  section  and  in  Table  2  are  not  the  only  ways  to  identify  the  parameters.  The  parameters  may  be 


Figure  14:  Relationship  between  the  stress-strain  response  of  SMA  and  the  initial  estimate  for  the  relative 
stress  g r  and  final  switching  stress  gr. 
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Strain 


Figure  15:  Relationship  between  the  stress-strain  response  of  SMA  and  the  initial  estimate  for  the  relative 
stress  or. 


estimated  using  a  variety  of  tests;  however,  the  summary  presented  in  Table  2  attempts  to  minimize  the 
experiments  necessary  to  identify  the  parameters.  Using  these  techniques  requires  five  experiments,  listed  in 
Table  3.  If  a  resistance  model  is  unnecessary,  then  only  four  experiments  are  required. 

4  Parameter  Optimization 

Two  different  optimization  cases  are  presented  and  compared.  In  one  case,  only  the  densities  vr(<jr)  and 
vi(oi)  are  optimized  while  the  model  parameters  are  considered  constants  identified  from  the  experimen¬ 
tal  data.  In  the  second  case,  the  model  parameters  are  considered  initial  estimates  and  included  in  the 
optimization  algorithm  with  the  densities.  Both  cases  uses  MATLAB’s  Isqnonlin  algorithm. 

4.1  Densities  Only 

Optimizing  the  density  functions  (11)  and  (12)  involves  finding  the  coefficients  OLk,m  and  Pk  that  minimize 
the  sum  of  squared  errors 


k= 1 


(34) 


k= 1 


between  the  measured  strain  ek  and  predicted  strain  e  (cr  (tk>)  ,  T  (tfc)),  where  the  design  variables  are 

q  =  [ai,**-  ,OLKaiMa,Pirmm  iPKp\T- 


The  analytical  Jacobian  is  derived  to  decrease  computational  time  for  the  optimization  algorithm,  as  nu¬ 
merical  approximations  of  the  Jacobian  require  additional  function  evaluations  and  decreases  accuracy.  The 
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Table  2:  Data-driven  techniques  to  determine  initial  parameter  estimates. 


Parameter 

Data 

Source 

Ea 

Stress-strain  data  above  Tas 

Slope  of  the  high-temperature  stress-strain  response 
at  small  strain  values;  slope  of  line  ab  in  Fig.  6. 

Em 

Stress-strain  data  above  Tas 

Slope  of  the  high-temperature  stress-strain  response 
at  high  strain  values;  slope  of  line  cd  in  Fig.  6. 

st 

Stress-strain  data  above  Tas 

Intersection  of  the  high-temperature  martensite 
load-line  with  the  zero-stress  axis;  point  e  in  Fig.  6. 

°L 

Stress-strain  data  above  Tas 

Stress  midpoint  of  a  high-temperature  major  loop; 
point  /  in  Fig.  7. 

tl 

Stress-strain  data  above  Tas 

Temperature  at  which  cfl  is  estimated. 

Act  j1 

Stress-strain  data  at  two  differ¬ 
ent  temperatures  above  Tas 

Difference  in  the  midpoint  of  the  major  loops  at  two 
different  temperatures;  Eqn.  (16). 

t 

Creep  data  above  Tas  at  a  fixed 
strain  near  St 

Eqn.  (22). 

V 

Creep  data  above  Tas  at  a  fixed 
strain  near  1% 

Eqn.  (23). 

b 

Stress-strain  data  above  Tas 

Determine  the  95%  confidence  interval  over  which 
switching  occurs  (Fig.  13)  and  set  equal  to  a  a  ±  2  b 
or  fit  an  estimate  for  a  normal  density  to  the  exper¬ 
imental  data  (Eqn.  (31)). 

P R 

Stress-strain  data  above  Tas 

Estimate  the  difference  in  loading  and  unloading  at 
the  midpoint  of  the  major  loop  (Eqn.  (29))  or  fit  a 
log-normal  density  to  the  experimental  data  using 
the  mean  and  variance  (Eqn.  (32)). 

C 

Stress-strain  data  above  Tas 

Determined  from  the  95%  confidence  interval  when 
vr(<Jr)  and  z/j(crj)  stop  interacting  on  unloading 
(Eqn.  (30))  or  fit  a  log-normal  density  to  the 
experimental  data  using  the  mean  and  variance 
(Eqn.  (33)). 

Pmo 

Resistance  at  a  low  temperature 
and  St 

Eqn.  (24) 

Pao 

Resistance  at  a  high  temperature 
and  zero  strain 

Eqn.  (26) 

(Am 

Resistance-temperature  data  at 
low  temperatures  and  zero  stress 
or  resistance-strain  data  at  two 
temperatures 

Fig.  12(a)  and  Eqn.  (27)  or  Fig.  12(b) 

OLA 

Resistance-temperature  data  at 
low  temperatures  and  zero  stress 
or  resistance-strain  data  at  two 
temperatures 

Fig.  12(a)  and  Eqn.  (28)  or  Fig.  12(b) 

V 

Resistance-strain  data  at  low 
temperature 

Solve  Eqn.  (25) 
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Table  3:  Ideal  experiments  for  identifying  model  parameters. 


Description 

Stress-resistance-strain  data  at  temperature  T\  >  Tas  with  full  major  loop. 
Stress-resistance-strain  data  at  temperature  X2  >  T\  with  full  major  loop. 

Creep  data  at  temperature  T  >  Tas  and  6  ~  et-  quickly  load  the  wire  to  approximately  St  and 
then  hold  constant  while  recording  stress. 

Creep  data  at  temperature  T  >  Tas  and  e  ~  1%. 

Resistance-strain  data  at  temperature  X3  <  Tas  :  heat  the  wire  under  zero  stress  and  then  cool; 
load  to  fully  M+  and  then  unload. 


Jacobian  corresponding  to  (34)  is  given  by 


de(t 1) 

5e(t1) 

de{tx) 

de(P) 

da  1,1 

daKa,Mo> 

dp! 

d0Kp 

de(tNk) 

de(tNk ) 

de(tNk ) 

de(tNk) 

da  1,1 

daKa,Ma 

dp! 

d0Kfj  _ 

where 

~  (fk\  Ni  Ni  1 

-k—  =  EE  (  _  (^i)  ui  (aij) 1  (a(tk)  +  °ij>T(tk)><rRi)  wiwo  +  e  p(ifc),T(ifc))  ) 

uam  i=1  j=1  v  2^1  al  7 

r\  fj-k\  Ni  Nj  . 

iC=EE(  -  vR  {aRi)ipm  (aij)  £  ( a(tk )  +  <jIvT{tk),VRi)  u>iWj  +  ^  T(tk))  J . 


(35) 


(36) 


4.2  Densities  and  Model  Parameters 

Optimizing  the  densities  and  model  parameters  involves  two  minimization  routines.  In  the  first  step,  the 
sum  of  squared  errors 


1  Nk.  1  Nk.  /  2 

p  m  =  2  E (C2  =  2  E  (ik  - £  )  (37) 

z  k=  1  Z  k= 1 

between  the  measured  strain  and  predicted  strain  is  minimized.  Here  the  design  variables  are 

Qe  =  [oq,i,  •  •  •  ,  oiKa,MaiPii  •  ‘  •  5  Pkp,Ea,  Em,cjl ,  AaT,  t ,  V]T . 

Note  that  we  do  not  include  the  transition  temperature  Tl  due  to  its  coupling  with  <jl.  The  determination 
of  only  one  of  the  two  is  necessary. 

In  the  second  step,  the  sum  of  squared  errors 

pm  =  2  Ee^  kf  =  2  E  0  (*(**)>  r(**)) )  (38) 

z  fe=i  z  fe=i 

between  the  measured  resistance  Clk  and  the  predicted  resistance  Q  (a(tk),  T(tk))  is  minimized.  The  design 
variables  in  this  step  are 

On  =  [poA,PoM,MA,aM,v]T. 
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Again,  the  analytical  Jacobian  is  derived  to  decrease  computation  time.  The  partial  derivatives  corre¬ 
sponding  to  (37)  and  (38)  are  given 


de£{tk) 

dcXm 

de£(tk) 


:  See  Equation  (35) 
:  See  Equation  (36) 


d/3m 

de£(tk)  xA(tk)a(tk) 


dEA 


dEM 

de£(tk)  dxM+  ( <x{tk)  <J  (tk) 


de£(tk)  _  xM+(tk)  +xM-(tk) 


E 2 


Octl  9(Jl  V  Em  Ea 

de£(tfc)  _  _&cm±  ( °El  _  a(tk) 
dA cft  3 Act t  \  Em  Ea 

de£(tk) 

— ^ -  =  2A 

OEt 

de£{tk)  dxM+  ( cr(tk)  a(tk) 


+  £t 


+  £t 


dxM-  ( &(tk) 


drx  dr  \  EM  EA 

de£(tk)  dxM+  ( cr{tk)  a(tk) 


+  £r  }  — 


dV 


dV  V  Em 

ffe 


Ea 


■£t  I  ~ 


daL 

dxM-  , 

\  Em 
(a(tk) 

1  1 

dAcr  x 

V  Em 

Ea 

' M-{tk ). 

dxM-  , 

(a(tk) 

a(tk) 

dr 

\  Em 

Ea 

dxM-  - 

U(tk) 

a(tk) 

dV  1 

V  Em 

Ea 

■  £t 


■  St 


—  St 


—  St 


<9en(tfc)  L(  £e)  r,fcUl  /\T(tk\) 

dp0A  A(s(tk))  )  l1  +  aA^T{t  )) 


deq(tfc) 

dpoM 


L(se) 

'  A  (s(tk)) 

deq(tfc) 

daA 


(■ XM+(tk )  +  XM-(tk))  (l  +  «MAT(tfe)) 

L(£e) 


A{e{tk)) 


xA(tk)p0AAT 


=  ~JWY)  (W**)  +*M-(**))  P0MAT(tk) 
den(tk)  _  2L(se)  psrna(tk)ee{tk) 


du 


ttrI  (i  -  vse{tk)y 


The  partial  derivatives  of  the  phase  fractions  with  respect  to  ox,  A<tt,  t,  and  V  are  found  by 
ating  the  phase  fraction  evolution  ODEs  (4)  and  switching  the  order  of  the  derivative,  yielding 


=  -?p±AXM+{t) -p+a^E  +  §g±XA{t)_PA+  (^±0)  + 

=  _§p=AXM_{t)  _p+A^it)  +  ^XA{t)_pA_  + 


(39) 

(40) 

(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 
different  i- 

(53) 


where  7  =  ox,  A<tt,  t,  and  V.  The  derivatives  of  the  transition  rates  are 


dpap 

dr 


1 

Pa  [3  , 

T 
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and 


dpa0 

AGaf3 

dV 

kbT(t)PaP ’ 

9pap  

V  dAGa0 

doL 

kbT(t)Pa0  do  L  ' 

dpa/3 

V  dAGa0 

d/S.(7x 

kbT(t)Pai3  dAoT 

The  Gibbs  energy  derivatives  and  are  found  by  differentiating  (1).  Finally,  the  partial  deriva¬ 

tives  of  the  phase  fractions  are  found  by  integrating  (53). 

5  Experimental  Setup 

As  shown  in  Figure  16,  the  SMA  experiments  were  conducted  in  a  temperature-controlled  water  bath  to 
ensure  constant  temperature  conditions.  The  water  temperature  was  regulated  using  a  heating  element.  A 
low-temperature  FLEXINOL  wire  ( Tas  =  343  IT,  diameter  =  0.005  in.  (0.127  mm),  Dynalloy  Inc.,  Tustin, 
CA)  was  used.  The  displacement  of  one  end  of  the  wire  was  controlled  with  a  Micro-Drives  motor  and  BOP 
20-10m  operational  power  supply  (Kepco  Inc.,  Flushing,  NY)  and  measured  with  a  mDMS-D170  infrared 
displacement  sensor  (Philtec  Inc.,  Annapolis,  MD).  For  resistance  measurements,  leads  were  connected  to  the 
wire  and  a  small  current  (50  mA)  was  applied  using  a  6454A  operational  power  supply  (Agilent  Technologies, 
Santa  Clara,  CA).  The  voltage  was  monitored  in  order  to  determine  the  resistance,  and  the  force  in  the  SMA 
wire  was  measured  using  a  MLP-10  load  cell  (Transducer  Techniques,  Temecula,  CA). 

The  experiments  were  conducted  over  a  range  of  temperatures:  308,  318,  328,  338,  343,  348,  353,  and 
358  K.  The  zero-strain  lengths  in  the  wire  were  determined  from  its  austenitic  length.  At  temperatures  above 
343  K  (the  austenitic  start  temperature),  the  zero-strain  length  corresponds  to  the  point  where  a  force  in  the 
wire  is  initially  detected.  The  length  at  343  K  is  used  for  all  tests  below  this  temperature,  as  the  wire  can 
exist  in  a  blend  of  all  three  phases  below  the  austenitic  start  temperature  due  to  material  inhomogeneities. 
A  specific  procedure  described  below  was  followed  for  each  test  to  ensure  a  known  initial  condition. 

1.  Heat  the  wire  to  358  K  and  strain  it  to  3.5%  strain.  Apply  0.6  A  current  and  measure  the  force  (To). 

2.  Set  the  desired  water  temperature. 


Figure  16:  Experimental  setup. 
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3.  Displace  the  wire  to  3.5%.  Increase  the  current  in  the  wire  using  the  Agilent  power  supply  until  the 
force  equals  the  force  achieved  in  the  358  K  wire  at  3.5%  and  0.6  A  (Fq).  Maintain  this  current  while 
the  displacement  is  adjusted  to  the  zero- strain  length. 

4.  Set  the  current  to  zero  and  let  the  wire  cool  for  3  minutes  while  the  water  is  stirred  to  reach  equilibrium. 

5.  Strain  the  wire  to  7.0%  and  back  to  0.0%. 

6.  Repeat  steps  2-5  three  times. 

7.  Go  to  the  next  temperature,  starting  with  step  2. 

6  Model  Validation  Results 

6.1  Estimation  of  Model  Parameters 

Initial  estimates  of  the  model  parameters  are  obtained  using  the  techniques  described  in  Section  3.  The 
experimental  stress- strain  data  at  two  different  temperatures  is  shown  in  Figure  17(a)  and  the  creep  data  at 
two  different  strain  values  is  shown  in  Figure  17(b).  The  resistance-strain  data  used  to  estimate  pmo  and  v 
is  shown  in  Figure  17(c),  and  the  data  used  to  estimate  pAo  is  shown  in  Figure  17(d).  The  resistance-strain 
data  at  two  temperatures  is  shown  in  Figure  17(e)  and  used  to  estimate  q^m-  An  initial  estimate  for  a  a 
could  not  be  obtained  due  to  a  decrease  in  resistance  as  temperature  increased,  as  shown  in  Figure  17(f). 
This  phenomenon  is  likely  due  to  more  of  the  wire  converting  to  austenite  as  the  temperature  increases  and 
the  fact  that  the  nominal  resistivity  of  austenite  is  lower  than  martensite.  Therefore,  the  estimate  for  olm  is 
used  for  a  a-  The  initial  estimates  for  all  the  model  parameters  are  listed  in  Table  6. 

Two  of  the  parameter  estimates  listed  in  Table  6  deviate  from  expected  and  published  values.  The 
austenitic  elastic  modulus  is  considerably  lower  than  other  published  values  (70  GPa).  The  lower  value  may 
be  due  to  compliance  in  the  experimental  setup  or  slipping  of  the  wire.  Additionally,  the  Poisson’s  ratio  is 
approximately  twice  the  value  expected  for  the  material.  This  may  be  due  to  error  in  measuring  resistance 
or  the  power  supply  leads  being  affected  during  loading. 

Initial  estimates  for  the  relative  and  interaction  stress  densities  are  found  by  determining  the  experimental 
mean  and  standard  deviations  for  the  relative  stress,  148.1  MPa  and  30.5  MPa,  respectively,  and  standard 
deviation  for  the  interaction  stress,  35.4  MPa.  The  relative  stress  parameters  are  then  converted  to  log¬ 
normal  form  using  (32)  and  (33).  This  method  is  used  because  it  was  found  that  it  produces  a  better  fit 
versus  estimating  the  parameters  using  the  95.4%  confidence  intervals  (Figures  13  and  14).  The  relative 
stress  and  interaction  stress  densities  at  353  K  are  shown  in  Figure  18. 

6.2  Densities  Only 

The  parameters  for  the  density  basis  functions  (13)  and  (14)  are  chosen  a  priori  using  the  initial  estimates  for 
the  relative  stress  mean  and  standard  deviation.  Three  different  values  for  the  mean  /i#m  are  used:  E  (<r^) 
and  E  (<j^)  db  75MPa.  Two  different  values  for  the  standard  deviation  Ck  are  used:  the  initial  estimate  (30.5 
MPa)  and  five  times  this  value  (152.5  MPa).  A  combination  of  these  values  produces  six  basis  functions  for 
the  relative  stress  density  vr{gr).  The  basis  functions  for  the  interaction  stress  density  are  parameterized 
with  five  different  values:  0.255,  0.55,  5,  25,  and  45. 

The  stress-strain  data  at  338  and  353  K  is  used  to  optimize  the  density  coefficients  and  /?£.  The 
experimental  data  at  338  K  is  used  instead  of  the  lowest  temperature  data  (308  K)  due  to  the  two-way  shape 
memory  effect  seen  in  the  FLEXINOL  wires  [5].  Since  the  model  does  not  currently  quantify  this  effect,  the 
experimental  data  immediately  below  the  austenitic  start  temperature  is  used. 

Comparisons  between  the  experimental  stress-strain  data  and  the  model  are  shown  in  Figure  19,  and 
comparisons  between  the  experimental  resistance-strain  data  and  the  model  are  shown  in  Figure  20.  The  fit 
model  (338  and  353  K)  is  shown  in  Figures  19(a)  and  (b).  All  other  results  (308,  328,  343,  and  348  K)  are 
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predicted  model  results ,  including  all  the  resistance-strain  plots  in  Figure  20.  The  final  cost  (sum  of  squared 
error)  for  the  fit  model  is  0.0422  (or  a  mean  squared  error  of  3.01  x  10-5). 


(c) 
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Figure  17:  Identification  of  model  parameters:  (a)  stress-strain  data  at  348  K  and  358  K;  (b)  creep  data  at 
5.0%  and  2.0%  strain;  (c)  low-temperature  resistance-strain  data  for  estimating  pmo  (point  A)  and  v  (change 
over  A  to  B);  (d)  high-temperature  resistance-strain  data  for  estimating  pAo  (point  C);  (e)  resistance-strain 
data  at  two  temperatures  for  estimating  aMi  (f)  resistance-strain  data  at  two  temperatures  in  the  austenitic 
regime. 


(a) 


(b) 


Figure  18:  Identification  of  density  parameters:  (a)  distribution  of  experimental  relative  stress;  (b)  distribu¬ 
tion  of  experimental  interaction  stress. 
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(d)  (e)  (f) 

Figure  19:  Stress  versus  strain  results  for  optimized  densities  and  model  parameters:  Comparison  between 
model  fit  and  experimental  data  at  (a)  338  K  and  (b)  353  K;  Comparison  between  model  prediction  and 
experimental  data  at  (c)  308  K,  (d)  328  K,  (e)  343  K,  and  (f)  348  K. 


Figures  19  and  20  reveal  a  number  of  interesting  results.  As  expected,  the  model  fails  to  accurately 
quantify  the  lowest  temperature  (308  K)  response  shown  in  Figure  19(c),  where  the  material  is  undergoing 
a  transition  from  a  blend  of  M+/M—  to  M+.  The  model  predicts  a  much  higher  transition  stress  than 
observed  in  experimental  data.  While  the  data  is  not  used  in  the  optimization  algorithm,  it  is  worth  noting, 
that  the  lowest  temperature  data  has  the  highest  mean  squared  error  (MSE):  2.72  x  10-4.  The  other  predicted 
model  results  (328,  343,  and  348  K)  have  a  MSE  on  the  order  of  the  fit  model:  5.23  x  10— 5 ,  2.24  x  10-5,  and 
2.71  x  10  5,  respectively.  These  results  indicate  the  model  is  accurately  quantifying  the  material  behavior 
without  over-fitting  the  data. 

The  resistance-strain  data  in  Figure  20  shows  that  the  model  matches  the  experimental  data  best  at 
high  temperatures;  see  Figures  20(b),  (c),  and  (f).  As  actuators  operate  at  higher  temperatures  (in  order  to 
recover  strain  above  Tas )  the  results  indicate  that  the  model  may  be  suitable  for  self-sensing  applications, 
where  the  resistance  in  the  wire  is  used  to  predict  strain  [7,8].  However,  operating  in  low  temperatures  may 
be  problematic  for  these  applications,  as  multiple  strain  values  correspond  to  a  single  resistance  measurement 
(strains  below  2%  in  Figures  20(c)  and  (d)). 

The  optimal  density  coefficients  are  listed  in  Tables  4  and  5,  and  the  corresponding  densities  are  shown 
in  Figure  21.  The  relative  stress  density  is  shown  in  Figure  21(a),  and  the  interaction  stress  density  is  shown 
in  Figure  21(b).  The  optimal  density  coefficients  reveal  that  the  poly  crystalline  behavior  is  largely  modeled 
using  the  interaction  stress.  The  highest  coefficient  is  on  the  largest  standard  deviation  (141.60  MPa).  The 
model  averages  the  transition  behavior  without  exactly  quantifying  the  transition  stresses.  As  shown  in  the 
high  temperature  stress-strain  results  (Figures  19(b),  (e),  and  (f)),  the  model  matches  the  midpoints  of  the 
loading  and  unloading  curves,  but  fails  to  quantify  the  initial  and  final  transition  points  (around  1%  and 
4%  strain  in  the  figures).  This  discrepancy  is  likely  due  to  errors  in  the  initial  estimates  of  the  other  model 
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Figure  20:  Resistance  versus  strain  results  for  optimized  densities:  Comparison  between  model  fit  and 
experimental  data  at  (a)  338  K  and  (b)  353  K;  Comparison  between  model  prediction  and  experimental  data 
at  (c)  308  K,  (d)  328  K,  (e)  343  K,  and  (f)  348  K. 


ot  (MPa] 

(b) 


Figure  21:  Optimized  densities:  (a)  relative  stress  density  vr(cfr );  (b)  interaction  stress  density 


parameters,  which  motivates  the  need  for  optimizing  these  estimates. 

The  kernel  number  does  not  affect  the  computational  efficiency  of  the  model,  since  the  densities  and 
quadrature  weights  can  be  pre-computed  once  the  optimal  densities  are  known;  however,  the  optimal  coeffi¬ 
cients  listed  in  Tables  4  and  5  reveal  that  some  of  the  kernels  are  unnecessary.  The  kernels  with  coefficients 
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Table  4:  Optimized  coefficients  and  parameters  for  vr(<Jr)  for  optimization  of  the  densities  only. 


m 

k 

& k,m 

M Rm 

Ck 

1 

1 

0.000 

18.03 

0.159 

1 

2 

0.004 

17.80 

0.621 

2 

1 

0.324 

18.79 

0.041 

2 

2 

0.035 

18.72 

0.191 

3 

1 

0.000 

19.22 

0.019 

3 

2 

0.636 

19.18 

0.089 

Table  5:  Optimized  coefficients  and  parameters  for  z//(cr/)  for  optimization  of  the  densities  only. 


k 

Pk 

bk  (MPa) 

1 

0.000 

8.85 

2 

0.000 

17.70 

3 

0.329 

35.40 

4 

0.002 

70.80 

5 

0.669 

141.60 

less  than  0.005  can  be  neglected.  Alternatively,  the  kernels  with  large  optimal  coefficients  can  be  perturbed 
to  create  new  kernels,  and  the  optimization  algorithm  can  be  run  again. 

6.3  Densities  and  Model  Parameters 

For  the  optimization  of  the  densities  and  model  parameters,  the  same  density  basis  functions  are  used;  the 
values  for  (ir^  and  Ck  listed  in  Table  4  and  the  bk  values  listed  in  Table  5.  Since  the  strain  model  and 
resistance  model  parameters  are  optimized,  both  the  stress-strain  data  and  resistance-strain  data  at  338  K 
and  353  K  are  used.  The  stress-strain  model  parameters  and  densities  are  optimized  first.  Using  the  optimal 
values  for  the  strain  model,  the  resistance  model  parameters  are  optimized  next. 

Comparisons  between  the  model  and  experimental  data  are  shown  in  Figures  22  and  23.  The  fit  model 
(338  and  353  K)  is  shown  in  Figures  22(a)  and  (b)  and  Figures  23(a)  and  (b).  The  other  figures  compare  the 
predicted  model  to  the  experimental  data.  The  final  sum  of  squared  errors  for  the  stress-strain  optimization 
is  0.0060  (or  a  MSE  of  4.29  x  10  6).  The  sum  of  squared  errors  for  the  optimization  of  the  densities  and  model 
parameters  is  85%  lower  than  optimization  of  the  densities  only.  A  comparison  of  Figures  19(a)  and  (b)  and 
Figures  22  (a)  and  (b)  reveals  that  optimizing  the  model  parameters  allows  the  model  to  accurately  capture 
the  initial  and  final  transition  regions  from  austenite  to  martensite,  instead  of  averaging  the  response. 

Even  optimizing  the  model  parameters  fails  to  quantify  the  low-temperature  stress-strain  response  shown 
in  Figure  22(d).  The  MSE  for  the  308  K  data  is  2.55  x  10— 4 ,  similar  to  the  results  for  optimizing  densities  only. 
The  predicted  model  MSEs  compare  favorably  to  the  optimization  of  densities  only:  9.53  x  10— 5 ,  1.56  x  10— 5 , 
and  2.22  x  10-5  for  328,  343,  and  343  K,  respectively.  Whereas  the  fit  model  MSE  is  significantly  lower  for  the 
optimization  of  densities  and  model  parameters,  the  predicted  model  MSEs  are  the  same  order  of  magnitude 
as  the  optimization  of  densities  only,  suggesting  that  the  initial  estimates  for  the  model  parameters  are 
accurate.  The  optimized  model  parameters  are  listed  in  Table  6  along  with  the  initial  estimates. 

The  optimal  density  coefficients  are  listed  in  Tables  7  and  8.  The  corresponding  densities  are  shown  in 
Figures  24(a)  and  (b).  Comparing  these  densities  to  the  previous  ones  (Figure  21)  indicates  that  the  standard 
deviations  for  the  overall  densities  are  lower,  especially  for  the  interaction  stress  density.  Optimizing  all  the 
parameters  reduces  the  reliance  on  the  interaction  stress  density  for  modeling  the  polycrystalline  behavior. 
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Figure  22:  Stress  versus  strain  results  for  optimized  densities  and  model  parameters:  Comparison  between 
model  fit  and  experimental  data  at  (a)  338  K  and  (b)  353  K;  Comparison  between  model  prediction  and 
experimental  data  at  (c)  308  K,  (d)  328  K,  (e)  343  K,  and  (f)  348  K. 


Table  6:  Comparison  between  initial  estimates  of  SMA  model  parameters  and  optimal  values. 


Variable 

Initial  estimate 

Optimal  value 

Units 

Ea 

35.0 

35.0 

GPa 

Em 

28.0 

27.2 

GPa 

275 

295 

MPa 

AcrT 

10 

9.2 

MPa/K 

£T 

4.6 

4.44 

% 

T 

1.09 

1.78 

s 

V 

8.0xl0-26 

9.5x  10-26 

m 

POA 

8.8  x  10-7 

8.9xl0-7 

Ohm-m 

POM 

10.7  x  10-7 

10.7xl0-7 

Ohm-m 

OLA 

0.0003 

0.0001 

Ohm-m/K 

olm 

0.0003 

0.0001 

Ohm-m/K 

V 

0.67 

0.8 

- 

It  is  worth  noting  that  even  optimizing  the  model  on  the  308  K  data  fails  to  quantify  the  behavior 
observed  experimentally.  Using  the  low-temperature  data  in  the  optimization  algorithm  produces  a  MSE  of 
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Figure  23:  Resistance  versus  strain  results  for  optimized  densities  and  model  parameters:  Comparison 
between  model  fit  and  experimental  data  at  (a)  338  K  and  (b)  353  K;  Comparison  between  model  prediction 
and  experimental  data  at  (c)  308  K,  (d)  328  K,  (e)  343  K,  and  (f)  348  K. 


(a) 


(b) 


Figure  24:  Optimized  densities:  (a)  relative  stress  density  v\ (ctr);  (b)  interaction  stress  density  ^2(07)  • 


1.62  x  10-4.  This  value  is  lower  than  the  308  K  MSE  when  the  optimization  uses  the  338  and  353  K  data, 
but  it  is  still  an  order  of  magnitude  higher  than  the  high-temperature  data.  Additionally,  the  MSEs  for  the 
other  temperatures  increase. 
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Table  7:  Optimized  coefficients  and  parameters  for  vr(<jr)  for  optimization  of  the  densities  and  model 
parameters. 


m 

k 

(%k,m 

Ck 

1 

1 

0.010 

18.03 

0.159 

1 

2 

0.019 

17.80 

0.621 

2 

1 

0.306 

18.79 

0.041 

2 

2 

0.053 

18.72 

0.191 

3 

1 

0.271 

19.22 

0.019 

3 

2 

0.340 

19.18 

0.089 

Table  8:  Optimized  coefficients  and  parameters  for  z/j(crj)  for  optimization  of  the  densities  and  model 
parameters. 


k 

Pk 

bk  [MPa] 

1 

0.004 

8.85 

2 

0.031 

17.70 

3 

0.249 

35.40 

4 

0.075 

70.80 

5 

0.641 

141.60 

6.4  Minor  Loop  Comparison 

Minor  loop  data  was  obtained  on  the  experimental  setup  at  353  K  on  a  different  SMA  wire.  The  data  is 
collected  by  reversing  the  direction  of  stress  within  the  transition  region  identified  in  Figure  13,  either  on 
loading  or  unloading.  Example  experimental  data  that  quantifies  a  minor  loop  during  both  loading  and 
unloading  is  shown  in  Figure  25. 

Since  the  minor  loop  data  was  obtained  on  a  different  wire,  the  strain  model  parameters  (including  the 
densities)  are  optimized  again,  resulting  in  the  optimal  parameters  listed  in  Table  9.  The  optimal  density 
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Figure  25:  Example  minor  loop  data:  (a)  stress  vs.  time;  (b)  strain  vs.  time. 
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coefficients  are  listed  in  Tables  10  and  11.  Comparisons  between  the  fit  model  and  experimental  data  are 
shown  in  Figures  26(a)  and  (b).  Comparisons  between  the  model  prediction  and  experimental  data  are 
shown  in  Figures  26(c)  and  (d)  for  smaller  minor  loops. 

As  shown  in  Figure  26,  the  model  is  capable  of  quantifying  hysteresis  in  the  minor  loops.  However, 
the  model  predicts  greater  hysteresis  in  the  minor  loops  on  loading  (Figures  26(a)  and  (c))  than  unloading 
(Figures  26(b)  and  (d)).  This  effect  may  be  due  to  a  low  variance  in  the  interaction  stress  density.  Within 
these  minor  loops,  the  stress  does  not  reach  a  high  enough  value  to  begin  switching  kernels  back  to  M+. 

7  Conclusion 

In  this  paper,  we  have  presented  data-driven  techniques  to  estimate  the  parameters  for  the  homogenized 
energy  model  for  SMA.  The  model  is  described,  and  the  relationships  between  the  model  parameters  and  the 
material’s  response  are  discussed.  Techniques  are  developed  to  estimate  the  parameters  using  relatively  few 
experimental  tests.  These  estimates  are  used  in  two  manners.  In  the  first  method,  only  the  density  functions 
in  the  HEM  are  optimized.  In  the  second  method,  the  estimates  are  included  with  the  density  function 
coefficients  in  the  optimization  algorithm.  Optimizing  the  model  parameters  greatly  reduces  the  fit  model 
error;  however,  the  predicted  model  error  is  only  slightly  lower,  validating  the  techniques  for  identifying  the 
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Figure  26:  Stress  versus  strain  results  for  minor  loops:  Comparison  between  model  fit  and  experimental  data 
for  minor  loop  on  (a)  loading  and  (b)  unloading;  Comparison  between  model  prediction  and  experimental 
data  for  minor  loop  on  (c)  loading  and  (d)  unloading. 
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Table  9:  Optimal  model  parameters  for  minor  loop  data. 


Variable 

Optimal  value 

Units 

Ea 

35.0 

GPa 

Em 

28.9 

GPa 

265 

MPa 

/S.CFr 

8.8 

MPa/K 

4.44 

% 

T 

1.606 

s 

V 

4.2xl0-25 

m 

Table  10:  Optimized  coefficients  and  parameters  for  vr(<jr)  for  minor  loop  data. 


m 

k 

&k,m 

^Rrn 

Ck 

1 

1 

0.050 

18.03 

0.159 

1 

2 

0.014 

17.80 

0.621 

2 

1 

0.012 

18.79 

0.041 

2 

2 

0.004 

18.72 

0.191 

3 

1 

0.757 

19.22 

0.019 

3 

2 

0.164 

19.18 

0.089 

Table  11:  Optimized  coefficients  and  parameters  for  i for  minor  loop  data. 


k 

Pk 

bk  [MPa] 

1 

0.423 

8.85 

2 

0.064 

17.70 

3 

0.070 

35.40 

4 

0.138 

70.80 

5 

0.305 

141.60 

model  parameters.  Both  approaches  accurately  quantify  the  macroscopic  behavior  of  the  material  at  higher 
temperatures  (>  32 8K)  but  fails  to  quantify  the  response  at  low  temperatures  (30 8K).  Future  work  will 
investigate  changes  to  the  mesoscopic  model  to  quantify  the  effects  seen  at  low  temperatures. 

The  techniques  presented  here  relate  model  parameters  to  the  material’s  physical  response,  facilitating 
model  adoption  by  reducing  the  model’s  complexity.  Furthermore,  the  ideal  experiments  for  estimating  the 
parameters  are  presented.  Using  relatively  few  experiments  (five),  the  HEM  parameters  can  be  related  to 
the  experimental  results. 

In  addition  to  accurately  quantifying  the  macroscopic  actuator  behavior,  the  homogenized  energy  model’s 
computational  efficiency  makes  it  ideally  suited  for  control  algorithms  and  design  optimization.  The  model 
also  accurately  represents  the  resistance-strain  response,  providing  for  its  use  in  self-sensing  applications. 
Future  work  will  incorporate  the  model  into  real-time  control  and  sensing  applications. 
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